Macroscopic equations for pattern formation 
in mixtures of microtubules and motors 
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Inspired by patterns observed in mixtures of microtubules and molecular motors, we propose 
continuum equations for the evolution of motor density, and microtubule orientation. The chief in- 
gredients are the transport of motors along tubules, and the alignment of tubules in the process. The 
macroscopic equations lead to aster and vortex patterns in qualitative agreement with experiments. 
While the early stages of evolution of tubules are similar to coarsening of spins following a quench, 
the rearrangement of motors leads to arrested coarsening at low densities. Even in one dimension, 
the equations exhibit a variety of interesting behaviors, such as symmetry breaking, moving fronts, 
and motor localization. 

PACS numbers; 87.10+e, 87.15-v, 05.10-a 



I. INTRODUCTION 

The emergence of complex patterns and correlated 
fluctuations is a characteristic of many out of equilib- 
rium situations (l]J^. Living organisms provide many 
examples, ranging from the flocking of birds the 
social organization of bacterial colonies ||5|-[7| , to internal 
reorganization of a cell during division 

At the molecular level, motors and microtubules are 
frequently the ingredients responsible for construction 
and motion. Molecular motors are the proteins that con- 
vert chemical energy to mechanical energy, and have been 
extensively studied JlO|-|l^. Microtubules, consisting of 
a subunit protein called tubulin provide the scaffolding 
for many cell constructs, as well as "railways" for trans- 
port of proteins p^JT^ . In particular, microtubules have 
a polarity that provides a direction for the transport of 
motors. 

One of the many processes in which motors and mi- 
crotubules are involved is cell division. They are indeed 
the main ingredients of mitotic spindles which serve to 
move apart the duplicated chromosomes of eucaryotic 
cells 1^,^ . To elucidate some of the physical mechanisms 
involved, several in vitro experiments on mixtures of mi- 
crotubules and motors have been carried out [^-^. 
Even these simple mixtures result in interesting patterns: 
At an initial stage the microtubules form an aster with 
their "plus" ends pointing towards a center. The kinesin 
based beads move along the microtubules towards this 
center. In a conflned geometry, the aster pattern is then 
destabilized, giving way to a vortex in which the motors 
rotate around a center. At larger scales and in uncon- 
fined geometries, a variety of self-organized patterns are 
obtained upon varying the motor concentration. With 
increasing concentration, an array of vortices, a mixture 
of asters and vortices, a collection of asters, and bundles 
of microtubules emerge. 

The motivation of this paper is to describe some of 



the patterns observed in these experiments. There are 
in fact already several models of these phenomena in the 
literature, starting with the simulations reported in the 
original paper on the patterns |p^ . A two-dimensional 
model for the orientations of microtubules was used to 
produce inhomogeneous stripe patterns | [T^ . Another re- 
cent model introduces a convection-diffusion equation for 
motor density in the presence of the microtubule array. 
This model results in a density profile of motors in asters 



which decays as a power-law 19 1. 

Our approach is to take a macroscopic perspective, in- 
troducing two continuous fields m(r, t) and T(f, t) to de- 
scribe the local motor density and tubule orientation, 
respectively. The central input to these equations is that 
the motor density is transported along the tubule direc- 
tion, while the tubules are in turn aligned by the motors. 
In the in vitro experiments |l^], the physical origin of 
the latter is that the motor complexes can attach to two 
nearby tubules and their motion along the two provides 
a force that makes them parallel In writing such 

continuum equations for densities, we are following re- 
cent studies modeling the flocking of birds ||,^ , and the 
organization of growing bacteria (^-0]. Simulations of 
these equations reproduce asters and vortices in agree- 
ment with experiment. Analytical solutions also provide 
further insights on the patterns. For example, we find 
that the motor density is much larger close to the center 
of an aster, than in a vortex. The resulting increased 
strain energy on the tubules provides a driving force for 
asters to break off to form vortices, in qualitative agree- 
ment with experiments. 

The global evolution of the tubule pattern is sensitive 
to the initial density of motors. At high density, it is sim- 
ilar to coarsening of XY spins following a quench from 
high temperatures pl|-p5t . With periodic boundary con- 
ditions, the ultimate pattern is one of aligned tubules, 
with the motors going around in a uniform current. Such 
a pattern is not possible with closed boundary conditions, 



1 



which typically lead to a single vortex in the center. At 
lower densities, fluctuations play a strong role, and we 
observe a novel phenomenon of arrested coarsening, in 
which an inhomogeneous pattern of tubules freezes at 
some point in time. This occurs because the transport of 
motors produces regions in which the density of motors 
is very low (effectively zero). When such regions perco- 
late throughout the system, no further rearrangement of 
tubules is possible. 

To better understand the dynamics of coarsening and 
sizes of the arrested domains, we also simulated the equa- 
tions in one dimension. Even in this case, the equations 
exhibit a variety of interesting patterns that depend on 
the boundary conditions: (i) With periodic boundary 
conditions, there is a phase transition between a state 
with a uniform current of motors running along tubules 
aligned in one direction (at high motor densities); and 
one in which there is a localized cluster of motors moving 
at constant velocity around the system (at low densities). 
Note that both patterns correspond to a broken symme- 
try (of the two possible tubule orientations) in one di- 
mension. In contrast to the equilibrium Ising model, this 
symmetry breaking appears to persist in the presence 
of noise (mimicking finite temperatures), (ii) Reflecting 
boundary conditions lead to an oscillating front sweeping 
back and forth across the system; yet another solitonic 
solution to the equations, (iii) Closed boundary condi- 
tions give rise to initial coarsening, and eventual freezing 
of the tubules into domains. The domain size for frozen 
tubules depends on the value of the average motor den- 
sity. Unlike its two dimensional counterpart, the motor 
density continues to evolve after the tubules are frozen, 
and all motors are eventually localized in one cluster. 



II. MODEL 

We introduce the local motor density m{x, t), and the 
tubule orientation field T{x,t). The conservation of mo- 
tors leads to the continuity equation 



dm 



(1) 



where for the motor current, we shall assume the form 



D Vm + AmT. 



(2) 



Here, D is the diffusion constant for motors, while A is 
a coefficient describing their transport along the tubule 
direction T. These are the lowest order terms in an ex- 
pansion in m; higher order terms are expected and may 
become important at high motor densities. While tubules 
generally grow and fragment ||l^; their lengths can be 
stabilized by addition of taxol. We model this prefer- 
ence towards a particular length by an energy function, 
— -I- /3r^/4, similar to that used to describe spins 



in the Landau- Ginzburg equation. In the same spirit, we 

^ 2 

associate an energy cost of K{\/T) /2 with variations in 
the orientation of tubules. Extremizing such an energy 
function leads to 



dT 

'at 



aT- PT^T + V ■ {KVf). 



(3) 



Since the alignment of tubules is entirely by motors, we 
set K = 7m. Of course, since this is a non-equilibrium 
process, there is no a priori reason for the dynamics to 
originate from the minimization of an energy functional. 
A more systematic approach is to include all terms al- 
lowed by symmetries. In this case the equation for dT/dt 
may include terms jm'V'^T and 7'Vm • VT with different 
coefficients. For ease of description we shall start with the 
case 7 = 7', so that the dynamics is similar to minimizing 

an energy locally proportional to jm{x,t){\/T) /2. 

Finally, we can bring the equations into simpler form 
by rescaling to 



dm 
'dt 
df 
'dt 



V^m- V • (mf), 
Cf(l-T2)+V-(mVf), 



(4) 
(5) 



where we now measure length in units of y^(3/aD/A, 
time in units of /3D/ (aA^), motor density in units of 
-D/7, and the tubule density vector in units of \J Oij (5. 
The remaining parameter C is given by j3D/A^. 



III. SIMULATIONS 

We perform numerical simulations on a two- 
dimensional Lx L lattice, adapting the Crank-Nicholson 
scheme with the ADI operator splitting method [ p6| . 
The equations are discretized with spatial intervals of 
Aa; — Ay — 1, and time intervals of At — 10~^. At the 
edges of a finite system, we employ one of several possible 
boundary conditions: Reflecting boundaries have fixed in- 
ward pointing microtubules described by 



boundary 



(6) 



where n is the normal outward vector at the bound- 
ary. This discourages motors from approaching an edge. 
By contrast, with parallel boundary conditions the mi- 
crotubules are tangential to the boundaries, while closed 
boundary conditions place no restriction on the tubules. 
In all these cases, there is no current transporting motors 
outside the system. There is no restriction on the motor 
current when periodic boundary conditions are applied, 
although again the total number of motors is conserved. 

We start with an initial condition in which the mo- 
tor density is uniformly set to mo at all points, while 
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the tubule field has magnitude |To| = 10~^ and random 
orientations. After a transient period, the homogeneous 
configuration self organizes into patterns that depend on 
the value of average motor density mp, as well as the 
growth constant C. Figure |l| shows a mixture of vor- 
tices and asters which arises as the stationary pattern for 
mo = 0.01. Both asters (tubules pointing inward) and 
vortices (tubules going around) are clearly visible and 
randomly arranged throughout the system. The motor 
density now becomes inhomogeneous, with motors accu- 
mulating in the centers of vortices and asters [Fig. |^] . The 
asters are more visible and dominant as the initial motor 
density increases, as shown for mo = 0.15 in Fig. [H How- 
ever, higher densities lead to a single vortex as in Fig. ^ 
for mo = 0.5 I?!- 
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FIG. 1. A lattice of vortices and asters for mo=0.01 and 
C=100. Initially, the size of the tubule is 10"'^, and their 
directions are random. 
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FIG. 2. The profile of motor density corresponding to 
Fig. |l| The initial homogeneous density evolves to obtain 
peaks at the centers of asters and vortices. 
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FIG. 3. A lattice of asters for mo = 0.15 and C = 100. 
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FIG. 4. A large vortex for mo = 0.5 and C = 100. At 
sufficiently high densities of motors, one or several vortices 
are formed. 

The dynamics of formation of a large vortex from the 
random initial conditions is depicted in Fig. ^j. This fig- 
ure corresponds to reflecting boundary conditions, with 
C=10, and mo=0.15. Since the initial tubule length is 
smaller than unity, the first stages of evolution are the 
lengthening of tubules as depicted in Figs. ||(a) and||(b), 
for times t = 0.6 and t = I, respectively. During this 
stage the directions of the tubules do not change and 
remain randomly distributed. The next stage involves 
reorientations of the tubules. Since a uniform alignment 
is incompatible with the boundary conditions, an aster 
forms in the center as depicted in Fig. ||(c), for t ~ 120. 
Motors are now transported along the tubules and accu- 
mulate at the center of the aster. At longer times, the 
aster pattern gives way to a vortex as in Fig. ||(d) for 
t = 1200. The vortex pattern is stable, although its cen- 
ter may move around depending on the choice of bound- 
ary conditions ||2^ . This dynamics is consistent with 
the experimental observation that vortices form from the 
destabifization of asters IM . 
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FIG. 5. The evolution of a vortex from an aster at times 
t = 0.6 (a),t = l {h),t= 120 (c), and t = 1200 (d). The ran- 
dom initial pattern of tubules first grows in length, and then 
organizes as an aster. Motors accumulate at the center of the 
aster, and then circle around when it changes to a vortex. 



To model pattern formation in larger systems in which 
boundary effects are less important, we also performed 
simulations with closed boundary conditions. This 
change in the boundary condition does not qualitatively 
alter the pattern formation process. At low motor den- 
sities, we still observe a mixture of vortices and asters, 
followed by a collection of asters as mo is increased. How- 
ever, at large motor density, after formation of a large 
vortex at the center, motors also pile up at several points 
on the boundary, as indicated in Fig. ^. 

Equations (Q) and (^) are deterministic; the only 
stochasticity appearing through the initial conditions. 
However, randomness and noise are certainly present in 
the experimental situations. In particular, microtubules 
are known to constantly grow and shrink though a dy- 
namic instability ]l3| , p^, while asters are still observed 
under such conditions [p8[. To make sure that the pat- 
terns observed in our simulations survive the addition 
of noise, we also introduced a stochastic noise in the 
tubule evolution equation. We observed that the self- 
organized patterns are stable at small noise, but that 
sufficiently large noise causes a phase transition to ho- 
mogeneous mixtures (bundles of microtubules in a uni- 
form motor density). Alternatively, we observe that for 
a fixed amount of noise, patterns are destroyed at motor 
densities lower than a critical rric, but are qualitatively 
unchanged otherwise. (The value of nic is 0.005 if the 
noise distributed uniformly between —1 and 1.) 
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FIG. 6. The configuration of tubules for mo = 0.5, with 
closed boundary conditions (a) ; and the corresponding profile 
of motor density (b). Motors pile up at several points on the 
boundary, in addition to the interior. 



IV. ASTER AND VORTEX SOLUTIONS 



We can easily find analytical solutions to Eqs. and 
(^) that describe the motor density in the center of an 
aster or vortex. To this end, we look for stationary so- 
lutions, dt m = dt T = 0, with radial symmetry. In 
an aster the tubules are directed towards the center and 
T = — f, where f is the unit radial vector. Balancing 
the diffusive current —drm, with that transported along 
tubules gives an exponential form 

"T.astcr('') = m(0)exp(-r). (7) 

This exponential profile is indeed verified by the simula- 
tions as depicted in Fig. ^ 
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FIG. 7. The points represent the profile of motor density 
in an aster pattern for mo = 0.1 and C = 100. The dashed 
line is a least-square fit to the exponential decay 35.4 exp(— r). 

The tubules go around the center of a vortex, and 
T = 6, where 9 indicates the tangential unit vector. Mo- 
tors are then transported in a uniform circular current 
by the tubules. To ensure that there is no radial current 
of motors we need dr{rdrm) = 0, whose solution is 



TOvortox(7') = -Af log (r/i?) , 



(8) 



where i? is a long distance cut-off, of the order of the vor- 
tex size. A logarithmic fit to the simulated vortex motor 
density profile is shown in Fig. H. 



FIG. 8. The simulated profile of the motor density of 
a vortex for mo — 0.2 and C = 1. The dashed line is the 
least-square fit to the form, — 3.3 log(r/14.2). 

The profile of motors in quasi two dimensional asters 
has been studied theoretically and experimentally in 
Ref. jl^. In this study the tubules can grow to form 
quite long tracks, and their density also falls off away 
from the aster center. These differences in tubule behav- 
ior (as opposed to our case where tubules have uniform 
density and length) lead to predicted power law decays 
of the motor density profile. 

The aster and vortex configurations of tubules are re- 
lated to topological defects in the XY model. However, 
they are equivalent defects in the XY model as one can 
be deformed into the other through a 90° rotation. The 
presence and rearrangement of motors in our problem 



breaks this symmetry and the two configurations become 
inequivalent. In particular, the two defects have very 
different static energies E = J (Prm{r){VTY 12. In the 
aster, the motors are concentrated close to the center 
leading to a high strain energy. By contrast, the mo- 
tor density in a vortex in more uniform. Consequently, 
for the same number of motors, a large aster has much 
higher energy than a large vortex. Since the dynamics 
tends to minimize this energy, we have an explanation 
for why asters give way to the more stable vortices. Pre- 
sumably, finite size effects in smaller asters, of the order 
of the decay length implicit in Eq. (|^), account for their 
stability at small motor densities as in Fig. ^. 



V. ARRESTED COARSENING 

If the motor density is maintained at a uniform and 
fixed value, the dynamics of the tubules is identical to the 
coarsening of an XY system following a quench from high 
temperatures. This problem has been extensively stud- 
ied [pT|-^ and (up to logarithmic corrections) the typical 
length scale of the pattern coarsens as f ~ i^/^. In our 
case, the motors rearrange themselves in the landscape of 
the tubules and the coarsening scenario is modified when 
the motor density fluctuations become significant. 

At high motor densities there are enough motors left 
over after formation of asters and vortices to cause 
further rearrangements of the tubules, and coarsening 
continues towards the final pattern consistent with the 
boundary conditions. However, at low densities the mo- 
tors quickly migrate to the centers of asters and vortices. 
The little motor density left in the regions in between 
defects may then be too small to cause further realign- 
ment of tubules which became frozen. We call this phe- 
nomenon arrested coarsening of tubules. The limiting 
value of Too for the onset of arrested coarsening in fact 
depends on the growth constant of tubules C, as indi- 
cated by the "phase diagram" sketched in Fig. |^. 



aster-vortex y 


vortex 







FIG. 9. The phase boundary in the plane of mo and log C. 
To the right of the boundary the pattern coarsens to a large 
vortex. The arrested coarsening on the left side leads to frozen 
asters and vortices, with the aster pattern becoming more 
dominant as mo increases. 
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VI. ONE DIMENSION 



article, and will be presented elsewhere. 



To further understand the patterns, we examine the 
equations in one dimension, where more detailed simula- 
tions are possible. In particular, we consider scalar fields 
m{x,t) and T(x,t) evolving as 



dtm = dx^m - dx{Tm), 
dtT^C{T-T^)+d,{md,T). 



(9) 
(10) 



Such equations may be appropriate to describe the move- 
ment of myosin motors along quasi-one dimensional bun- 
dles of actin filaments which occurs in muscle contraction, 
and may be responsible for other types of cell motion 
Bundles of actin molecules can be formed in vitro 
in mixtures with other inert polymers. These bundles 
consist of random mixtures of actin molecules oriented 
in the two possible direction. However, in the presence 
of myosin motors and ATP, there is a sorting of polar- 
ity [^,|o| and active contraction of polar filaments . 
While the bundles are usually linear with open ends, then 
sometimes assemble into a ring It may also be pos- 
sible to artifically confine other mixtures of tubules and 
motors in quasi-one dimensional containers with various 
boundary conditions. In fact, we observe that the solu- 
tions to Eqs. and ( [lO| ) are quite sensitive to boundary 
conditions. Specifically, simulations show the following 
results: 

(i) Periodic boundary conditions correspond to plac- 
ing the system on a closed loop. We observe two types 
of symmetry breaking, depending on the initial density 
of motors. If toq is larger than a critical value of TOc 
(toc ~ 0.04 for C — 1), the tubule pattern coarsens until 
the tubules are all aligned in one direction (T = or 
— 1). This is accompanied by a uniform current of motors 
that goes around the system. This symmetry breaking 
is similar to that of an Ising model. However, unlike the 
equilibrium Ising model, the broken symmetry survives 
in the presence of random noise (simulating finite tem- 
peratures) in Eq. (10). Any domains of opposite spin 
formed due to randomness are annealed by a rush of mo- 
tors to the domain walls. The symmetry breaking in the 
presence of noise is due to advection term in Eq. (0) , and 
active transportation of motors along the tubules. Indeed 
we checked that when the advection term is removed, the 
presence of noise in Eq. ( p^ leads to finite domains. 

If mo is smaller than rUc, a novel final state emerges in 
which the motors gather together in a cluster that moves 
around the loop with a constant velocity. The tubules are 
again ordered in one direction, except near the cluster, 
where they briefly take the opposite alignment. Figure 
|To| shows that the configurations of the tubules and the 
profile of the motor density for mo = 0.02 < nic- We have 
verified by direct numerical integration that Eqs. (|^) and 
( p^ ) do indeed support such a solitonic solution. How- 
ever, the mathematical details are not in the spirit of this 




FIG. 10. The profiles of motor density (dashed line) and 
tubule orientation (solid line) in one dimension with periodic 
boundary conditions, for mo — 0.02 and C = 1. The whole 
pattern moves as indicated by the arrow, with a fixed speed. 

(ii) Reflecting boundary conditions were imposed by 
requiring the tubules at the edges to point inward, i.e. 
[r(0) = -1-1 and T{L) — —1]. There is an initial coarsen- 
ing period in which domains of -1-1 and —1 grow inward 
from the respective edges. However, in the final pattern 
the boundary between the -1-1 and —1 domains is not 
stationary, but sweeps back and forth across the system! 
The motors are again concentrated at the interface, with 
a solitonic profile as indicated in Fig. ^ 




FIG. 11. The profiles of motor density (dashed line) and 
tubule orientation (solid line) in one dimension with reflecting 
boundary conditions, for mo = 0.4 and C = 1. The cluster of 
motors moves as indicated by the arrow, but its direction is 
reversed at each boundary, as the profile oscillates back and 
forth. 

(Hi) Closed boundary conditions were also applied, 
with no restrictions on the value of T at the edges, but 
setting the outward motor current to zero. Above a crit- 
ical motor density, coarsening of tubules proceeds to a 
single domain of the size of the system, L. Following the 
tubules, the motors then pile up at one end of the system. 
The low density behavior in this case is similar to the ar- 
rested coarsening observed in the two dimensional case: 
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The initial growth of +1 and —1 domains is stopped at 
some point due to the local absence of motors necessary 
for continuing realignments. Figure |2| shows the average 
domain size, as a function of the average motor density. 
At the point when tubule evolution is stopped, the motor 
density has peaks at (+, — ) domain boundaries. How- 
ever, as time goes on, there is a slower ripening process 
in which the motors gradually diffuse against the unfa- 
vorable domains, and eventually aggregate at one point 
in the system, as in Fig. (It is not clear if this process 
is truly absent in two dimension, or merely takes too long 
to observe in simulations. There are similarities to the 
localization problem in which 2 is a critical dimension.) 
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FIG. 12. The average tubule domain size in a one dimen- 
sional system of length L with closed boundary conditions, as 
a function of average motor density mo for C = 1. 
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FIG. 13. The profiles of tubule orientation and motor den- 
sity (solid and dashed lines, respectively) for mo = 0.1 and 
C = 1, in a one dimensional system with closed boundaries. 

There are similarities in the above behaviors to pat- 
terns and phase transitions in other one dimensional non- 
equilibrium processes. A localization transition is ob- 
served in a population with diffusion, drift, and repro- 
duction in Refs. fS^-Q. This model of a bacterial colony 



living on an oasis in the presence of randomly fluctuat- 
ing wind, shows extended, localized, mixed states, as well 
as extinction, as a function of the average growth rate 
and convection velocity. In our model, the tubules (in 
their disordered state) provide a similar source of ran- 
dom advection. Ordered states are also observed in a 
lattice model with two species of particles in which the 
mobility of one species depends on the density of the 
other ll^,^. This model exhibits three phases: One 
with strong phase, a fluctuation-dominated phase, and 
another with uniform overall density. The couplings in 
our system between motors and tubules have similarities 
to the latter model, although the enforced conservation 
laws are different (there is no conservation of tubule ori- 
entations) . 

VII. DISCUSSION 

Given how little input is used to construct the macro- 
scopic equations (tubules transport motors and are 
aligned in the process), it is encouraging that many fea- 
tures of the experimental patterns are reproduced. As 
in the experiments, we observe arrays of asters and vor- 
tices, and large single vortices (formed from the break up 
of asters). However, at high densities of motors, the ex- 
periments lead to irregularly arranged bundles of tubules, 
a feature not present in our model. The absence of tubule 
bundles indicates the limitations of the macroscopic ap- 
proach. To reproduce the observed sequence of patterns, 
more physical input into the equations is necessary. For 
example, a potential modification is to include non- linear 
terms in the transported motor current to model the 
inter-particle interactions at high densities. 

In models of highway traffic , the current actually 
diminishes at high density, and can be approximated by 

JH = -Vm + f m e~™/™»='^ (11) 

where rTimax is a characteristic density for saturation of 
current. Since the alignment of tubules is intimately con- 
nected to the motor current, similar modifications should 
appear in Eq. (||). As discussed earlier, because of the 
non-equilibrium nature of the process, the equation for 
dtT does not have to be obtained from the minimiza- 
tion of an energy density. In fact, to describe saturation 
effects on alignment of tubules, we modified Eq. (||) to 

dtf = Cf (1 - T^) + Vm • Vf + me"™/'"-"''^ V^f + ry, 

(12) 

where is a stochastic noise. The sequence of tubule 
patterns obtained from these modified equations upon in- 
creasing initial motor density is depicted in Fig. We 
observe that the single vortex pattern becomes unstable 
at high densities, and we may regard the high density 
phase as representing bundles of tubules. We also find 
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that the patterns of vortices and asters are stable in the 
presence of noise in the following sense: If we average the 
whole configuration of tubules over time, the intermedi- 
ate regions between defects are washed away while the 
defects survive. 

The rich variety of behaviors observed in these simple 
equations are due to their noncquilibrium character. In 
the realm of equilibrium, for example, one dimensional 
systems at non-zero temperature are featureless and dis- 
ordered. Clearly non-equilibrium effects can lead to sym- 
metry breaking, and a rich interplay behaviors sensitive 
to boundaries. Biological systems can take advantage of 
such phenomena, and should indeed provide many inter- 
esting patterns in need of explanation. 
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(b) 



FIG. 14. The sequence of patterns obtained upon inclu- 
sion of saturation effects for mo = 0.05 (a), mo — 1 (b), and 
mo = 20 (c). 
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